Ion transport through confined ion channels in the presence of immobile charges 
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We study charge transport in an ionic solution in a confined nanoscale geometry in the presence of 
an externally applied electric field and immobile background charges. For a range of parameters, 
the ion current shows non-monotonic behavior as a function of the external ion concentration. For 
small applied electric field, the ion transport can be understood from simple analytic arguments, 
which are supported by Monte Carlo simulation. The results qualitatively explain measurements of 
ion current seen in a recent experiment on ion transport through a DNA-threaded nanopore (D. J. 
Bonthuis et. al, Phys. Rev. Lett, 97, 128104 (2006)). 
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Because of its central role in maintaining the home- 
ostasis of cells, ion transport throueh channels across cell 
membranes is of great importance ^ system 

with free ions, such as an aqueous solution, one might ex- 
pect the ion current / to increase with increasing external 
ion concentration c when a constant electric field is ap- 
plied. Surprisingly, in the presence of immobile charges 
fixed in the channel, the opposite may occur, with ion 
current decreasing with increasing c. For example, in the 
case of water-filled biological channels with strong ion 
binding sites, the ion conductance has been observed to 
reach a maximum and then decrease (or saturate) as c in- 
creases ; similar behavior is observed in DNA-threaded 
nanopores connecting two reservoirs 

An ion channel may be thought of as a thin hollow tube 
of length L where ions can enter or leave only through 
pores at the two ends. Because of the large difference 
in the dielectric constants of water (k^ w 80) and the 
membrane containing the channel {Km ~ 2), introducing 
an uncompensated ion into the channel requires overcom- 
ing an energy barrier due to the charge's self energy Us 
01 • The reason for this is that because ^ Km, an 
ion's electric field lines are concentrated inside the chan- 
nel over a length proportional to /ii/ K^/Km, where li is 
the shortest dimension of the channel [8|. The specific 
form Us takes depends on the nature of the channel. 
For a planar channel the electrostatic potential varies 
as U{r) ~ Inr for length scales l\\J k^ j Km > r > li, 
while for a linear channel U(r) ^ r. For channels which 
are relatively short and narrow, the larger dimension 
of the channel L ~ li^J K^j Km ^ h', this implies that 
the self-energies scale as \n{L/li)/li and L/lf in planar 
{li X L X L) and one-dimensional (/i x li x L) geome- 
tries respectively. For example, for a water-filled channel 
of dimensions Inm x Inm x bnm, Us is about 7 kT at 
T = 300K where k is the Boltzmann constant 

In addition to electrostatic interactions, one might in- 
quire as to the importance of hydrodynamics interac- 
tions. It is easy to see that hydrodynamic interactions 
are important only for systems which are much larger 
than some characteristic scale i?*. The length scale i?* 
may be estimated by comparing the electrostatic and the 



hydrodynamic forces between two ions separated by a 
distance r. The electrostatic interaction (in three dimen- 
■pj while the hydrodynamic force 

is Ih 



sions) is fp = -. ; — ^ 

lUdr^/r, 7 is the viscous drag coefficient, eg is 



the dielectric constant of the vacuum, tq is the radius 
of the ion, e is the charge of the ion and Ud is the ion 
drift velocity. Taking Ud = (E/j), where E is the elec- 
tric field acting on the ions in the channel, we obtain 
i?* — e/(47reoKuj7'oi?). For the experimental conditions 
of P, i?* w lOnm which is larger than the channel scale 
(the same result holds in two dimensions). In this paper 
our interest is in this regime, consequently we will ignore 
hydrodynamic interactions. 

Non-monotonic behavior in charged channels was pre- 
viously studied theoretically using a single vacancy model 
[lo| . under the assumption that the channel was strictly 
one-dimensional. A more recent study considered the 
ion current in a channel threaded with charged DNA, 
where the available space for ion motion was assumed 
to be effectively two-dimensional. In this case, the 
non-monotonic behavior was attributed to the two- 
dimensional specifics of the channel and the self-energy 
of the ions, and to a boundary layer effect at the edges 
of the channel 0]. 

In this paper we present a many-particle statistical 
model of interacting ions, and argue that, in the presence 
of fixed background charges inside the channel, the large 
self-energy of an individual ion is sufficient to give rise to 
a non-monotonic ion current / as a function of external 
ion concentration c. Our main result is that, irrespective 
of the effective channel dimension, there is a crossover 
temperature ~ Us/k, below which the ionic current 
may exhibit non-monotonicity. However, above T*, the 
current is a monotonically increasing function of c. Con- 
sequently, non-monotonic behavior can be observed only 
when Us is large enough for to be above the freezing 
temperature of water. For example, when U{r) l/r, as 
in large three-dimensional cavities, is much below the 
freezing temperature of water but when U{r) ~ Inr, the 
I vs c curve may have a minimum even at room temper- 
ature. In any case, for very high (or very low) density of 
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FIG. 1: Schematic phase diagram of ion transport in a channel 
of confined geometry at temperature T and volume V with 
Nb immobile background charges inside the channel. The 
thick line denotes the crossover temperature T, . The two in- 
sets are plots of ion current I versus external ion concentration 
c. 



background charges, / increases monotonically with c as 
is naively expected. This is summarized in Fig[TJ 

The above results are the consequence of two main 
competing mechanisms for ion transport: (1) Hopping 
current At low temperature, the fixed background 
charges are screened by counter- ions which thus re- 
side in close proximity to the background charges - one 
may think of the counter-ions as sitting 'on the sites' of 
the background charges. However, if one of the back- 
ground charges is not screened (a 'hole'), the screening 
counter-ion of an adjacent background charge can hop 
to it. Ih is approximately proportional to phipo ^ Ph), 
where ph is the density of holes, and po is some con- 
stant. Since ph decreases with increasing c, the ion cur- 
rent first increases, attains a maximum (at ph = po) and 
then decreases. (2) Bulk current /(,: Ions that are not 
strongly attached to any counterions will move more or 
less freely inside the channel, and, biased by the electric 
field, will contribute to the total current. /{, is a mono- 
tonically increasing function of c. The total ion current 
/ is sum of the hopping current Ih and the bulk current 
h, I = h + h. 

This intuitive picture for is supported by a simple 
model for driven diffusion, the partially asymmetric sim- 
ple exclusion process (PASEP) The PASEP consid- 
ers a one-dimensional lattice of sites, each of which may 
be either empty or occupied by a single particle. Particles 
may enter or leave the system at its ends, and a particle 
may hop to an adjacent site provided it is unoccupied. 
The parameters of the model are the rate of influx (a, 
7) and outflux (/?, S) of particles at the left and right 



ends, respectively, and the hopping rates between sites: 
g < 1 and 1, to the left and right, respectively (where the 
applied electric field may be thought of as the cause of 
asymmetry of the hopping rates). In the ion channel, a 
fixed charge screened by a counter-ion maps to an occu- 
pied site in the PASEP model, and an unscreened fixed 
charge maps to an unoccupied site in the PASEP. 

The phase diagram of the PASEP model has been fully 
elucidated (see for example, 12, [Si). If the incoming 
rates a and 7 are taken to be proportional to the outside 
concentration c, the behavior of current Ih can be imme- 
diately obtained using these results. It follows from 
that Ih ^ c for small c, and Ih ^ l/c for large c. At inter- 
mediate c, the current attains a maximum or a plateau. 
In the PASEP model, the various rates are taken to be 
constant, but, in reality, rates will depend on specific 
configurations of the system. Clearly the PASEP model 
cannot capture the appearance of the minimum in the / 
vs c curve. 

To understand this minimum we will consider a sta- 
tistical mechanical model of interacting ions in an ion 
channel where the channel is in contact with a reservoir 
of a fixed chemical potential /i and temperature T. For 
simplicity, we will consider a discrete model, where the 
positions of ions lie on a lattice. The kinetic energy of 
the ions is neglected, since ion motion in a fluid is over- 
damped. We assume that the electrostatic potential U (r) 
of a unit positive charge at position r inside the channel 
decays rapidly outside the channel. The Hamiltonian for 
a system of interacting charges of hardcore radius 7'o 



(1) 



where qi ~ ±1 is the charge of z-th ion, U{rij) is the 
interaction potential of ions i and j, whose separation 
is r,;j, N is the total number of ions and we denote 
Uo = U{rij = To); the self energy of an ion is given by 
Us = Uq/2. The definition of the Hamiltonian absorbs 
the chemical potential /i < 0, for simplicity assumed to 
be the same for both positive and negative charges, which 
is related to the fugacity z by z = exp{p/kT). Note that 
inserting a bound neutral pair (one one - charge) costs 
an energy —2p: because of cancellation there is no con- 
tribution from the first two terms in Eq. [1] The fugacity 
z controls the density of ions inside the channel. 

For small external electric fields, it is reasonable to 
assume that local thermal equilibrium is maintained. We 
thus include a constant external electric field Ex along 
the channel axis. Using J2z^j = - E,; if]^ 

Eq. [l]may be rewritten 

H = ^^q^qj[U{r^■i) -Uo] -E^q.,x, 

i^j i 

-{N+ -N^- NB fkT\n{zb) - NkT\n{z) (2) 
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where we have explicitly indicated the Nb fixed negative 
background charges, and where the sum is over all pairs 
of ions except those where both are background charges. 
Here and N_ are the total number of positive and 
negative mobile ions respectively, N = (N^+N^), Xj the 
X coordinate of i-th mobile ion, and zi, = exp[-~Us/kT]. 
For small Zf,, charge fluctuations in a finite channel are 
small, and (iV+ - ^_ - iVs) ~ [IH. 

We begin by considering the system at zero electric 
field; the charge distribution is then governed by the par- 
tition function Z = Y:^{l/N+\N_\)exjp{-Ho/T), where 
Hq is the Hamiltonian of Eq. [2] with E = Q and the sum 
is over all configurations. For small electric fields, the 
charge distribution will be essentially unchanged; we will 
use this to calculate the ionic current. 

Let us consider the fugacities z ~ zt, <C 1. Here we 
expand Z in powers of z and z^. Collecting leading order 
terms, we obtain 



hole density, is therefore given by 



Z 



Nb 



Zb + Z 

+e)(i/z^«+izb 



0(z 
-0{Vz 



Nb-2 4s 



Nb+2\ 



(3) 



where V is the channel volume measured in units of 
ionic volume. Numbering the terms on the right hand 
side of Eq. [31 wc may interpret them as follows: 
(1) one unscreened immobile charge (one hole). (2) all 
immobile charges are screened (no hole), (3) two un- 
screened immobile charges (two holes), (4) one excess 
positive or negative charge (apart from the screened 
backbone charges) and (5) one excess bound pair of pos- 
itive and negative charges. In terms (4) and (5), the 
factor V accounts for the possible placements of the ex- 
tra charges. Eq. [3] can be well-approximated by the 
first two terms alone for Zb z z^, where z, ~ 
■min{V~^^^,{zbNB/Vy^'^}- Now in this fugacity range, 
the probability Ph that there is exactly one hole can be 
written as Pji w NsZh/ {NsZh + z) using Eq. [3l 

Consider first the behavior at low fugacity, z < Zb, 
which we will call Region I. By examining Eq. [3] one can 
see that as long as z is not extremely small, the domi- 
nating configuration has one uncompensated background 
charge. The current flows by positive charges hopping 
from one background charge to another so that the hole 
moves from one end of the system to the other. The 
probability of having one hole in the system depends 
weakly on z in this regime. For the current to flow the 
hole must rccombinc with a charge from outside the pore, 
and this occurs with a recombination rate proportional 
to z/zf,. The current / is proportional to Ph time the 
recombination rate, giving I w Ih z. 

Now consider intermediate fugacities z^ < z z*. In 
this regime, which we term Region II, the probability Ph 
of having one hole goes as 1/z. Since for z > Zf, the 
recombination rate can be approximately taken as 1, the 
hole current, which in this regime is proportional to the 



(T X 



Zb 



{NBZb + z) 



(4) 



where cr is a constant related to the jump-rate of a hole 
from one site to another. In this fugacity range, there are 
no free bulk charges, so the total current I k, 1^ ^ 1/z, 
decreasing with increasing fugacity. 

Region III is the large fugacity limit z ^ z», where 
extra charges enter the system, although the background 
charges are already fully compensated. In this regime 
the current is clearly expected to increase with increasing 
fugacity. 

As a function of increasing z, we have the following: 
The current rises linearly in Region I, falls as 1/z in Re- 
gion II, and rises again in Region HI. Thus, it is the 
passage from Region I to II that determines the non- 
monotonic behavior. However, Region II may be unob- 
tainable - this happens when z* < z^. In this case. Region 
I crosses smoothly over to Region HI, and the ion current 
increases monotonically with z over its entire range. In 
other words, Region II is present only if T < T*, where 
= 2Us/k\n{V/NB)- Note that increases with the 
number of bound charges Nb ■ The above picture breaks 
down when the density of background charges is so high 
that ions can move freely (without hopping) from one 
background charge to another. This occurs when the 
typical distance between background charges is smaller 
than the screening length. Under such conditions we ex- 
pect the current to increase monotonically with fugacity. 

To support these simple arguments, we have per- 
formed Monte-Carlo simulations. For computational 
convenience, ions are only allowed to move in discrete 
steps on a square lattice. Ions can enter and leave the 
system only from two opposite surfaces, representing the 
pores of the channel. A site may accommodate at most 
one ion. We denote by AH the energy difference between 
configurations after and before a possible Monte-Carlo 
move, with H defined in Eq. [51 The simulation is car- 
ried out in the following way: At each time step, a lattice 
site is randomly chosen. If it is an empty boundary site, 
a positive (negative) charge is created with probability 
TOm{^,ie~~}. If the site is occupied, the charge is 
destroyed with probability min{\,e~^}. If the site is 
in the interior of the lattice and is occupied, its charge is 
moved to a randomly chosen unoccupied neighboring site 



with probability min{{ 



kT 



}. For _E = 0, the system 



eventually comes to equilibrium, while for E ^ 0, the sys- 
tem settles into a non-equilibrium steady state with a net 
ion current across the channel in the x-direction. 

Motivated by the experiment of Ref. which is ef- 
fectively two-dimensional, we performed a simulation on 
an L X L lattice using the above protocol, with the inter- 
action potential taken to be U{r) = {2e'^ / K.^^ro) \n{L / r) , 
where e is the electron charge, and tq ~ 0.35 nm [6[. 
An immobile linear array of equally spaced unit negative 
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T=150K 
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FIG. 2: Ion current I (in arbitrary unit) across a two- 
dimensional channel versus fugacity z is plotted for different 
temperatures with L = 30 x where diameter of A'^ 
ion being 0.26 nm, an electric field E = 5.4 meV/nm along 
X-axis, Nb = 6 negative immobile background charges. 
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FIG. 3: Current / (in arbitrary unit), scaled probability Ph 
of exactly one hole (scaling factor 0.00055) and scaled total 
number of ions/site n (scaling factor 0.00012) is plotted versus 
fugacity z with an electric field iJ = 1.1 meV/nm, tempera- 
ture r = 150K, L = 30 X djf where dx = 0.26nm is the 
diameter of a K ion, zt = 1.23 X 10"® and iVs = 6 negative 
immobile background charges. 



charges is placed on a line parallel to the x-axis in the 
middle of the channel, at y = L/2, to mimic the pres- 
ence of charged ss-DNA in the experiment. One should 
note that when L is large and Nb = 0, Eq. [2] is the 2D 
Coulomb gas Hamiltonian 14 1. 

The results of the simulations are presented in Fig. 
[5] where we plot the total ion current versus fugacity 
for different temperatures. As expected from the argu- 
ments presented above, the numerical results are qualita- 
tively different in two different temperature regimes. For 
T < T* « 300K, the ion current / first increases with z 
for small z, then reaches a maximum and subsequently 
decreases. Increasing z further, / reaches a minimum and 
then starts increasing with z. For T > the current / 
is a monotonically increasing function of fugacity z. The 
numerical value of T, given above is somewhat smaller 
than that given by 2Us / iklii{V/NB)) = 816A'. 

This is about what could be expected from such a simple 
argument. 

In Fig. [31 we plot numerical results for the ion current 
/, the probability Ph that the system has exactly one 
hole, and the average total number of ions per site n 
as a function of the fugacity where both and n are 
scaled suitably to relate to / for T < T^,. As can be 
seen at low fugacities {z < 0.01), single hole hopping is 
responsible for the ionic current. For larger fugacities, the 
number of free bulk charges increases, and the current, 
almost entirely due to flowing ions in the bulk, rises. For 
large self energies, as in the simulation, inserting a pair of 
positive and negative ions into the channel is much easier 
than inserting a single charge, so the number of unbound 



charges inside the system increases as 2^ for 2* < 2 <C 1. 
This is seen in Fig. [3l which shows a concomitant rise in 
ion current. 

Finally, it is worth noting the influence of the effective 
dimension of the system, which manifests itself in the 
functional form of the Coulomb interaction. In Monte- 
Carlo simulations of the same geometry (a 30 x 30 lat- 
tice with a linear array of negative immobile charges 
in the middle), but employing a Coulomb interaction 
U{r) = 1 /r, we found that the minimum in the current- 
fugacity plot can, in principle, also occur, but only at 
very low temperatures, of order T ^ 36 K which is clearly 
experimentally irrelevant. 
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